Draft version January 20, 2013 

Preprint typeset using 1^1^^ style emulatcapj v. 08/13/06 



O 



in 

6 



(N 
> 

p 

o 



X 



THE PROTOSTELLAR LUMINOSITY FUNCTION 
Stella S. R. Offner 

Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138 

AND 

Christopher F. McKee 

Physics Department and Astronomy Department, University of California, Berkeley, CA 94720 

Draft version January 20, 2013 

ABSTRACT 

The protostellar luminosity function (PLF) is the present-day luminosity function of the protostars in 
a region of star formation. It is determined using the protostellar mass function (PMF) in combination 
with a stellar evolutionary model that provides the luminosity as a function of instantaneous and final 
stellar mass. As in McKee & Offner (2010), we consider three main accretion models: the Isothermal 
Sphere model, the Turbulent Core model, and an approximation of the Competitive Accretion model. 
We also consider the effect of an accretion rate that tapers off linearly in time and an accelerating star 
formation rate. For each model, we characterize the luminosity distribution using the mean, median, 
maximum, ratio of the median to the mean, standard deviation of the logarithm of the luminosity, 
and the fraction of very low luminosity objects. We compare the models with bolometric luminosities 
observed in local star forming regions and find that models with an approximately constant accretion 
time, such as the Turbulent Core and Competitive Accretion models, appear to agree better with 
observation than those with a constant accretion rate, such as the Isothermal Sphere model. We 
show that observations of the mean protostellar luminosity in these nearby regions of low-mass star 
formation suggest a mean star formation time of 0.3±0.1 Myr. Such a timescale, together with 
some accretion that occurs non-radiatively and some that occurs in high-accretion, episodic bursts, 
resolves the classical "luminosity problem" in low-mass star formation, in which observed protostellar 
luminosities are significantly less than predicted. An accelerating star formation rate is one possible 
way of reconciling the observed star formation time and mean luminosity. Future observations will 
place tighter constraints on the observed luminosities, star formation time, and episodic accretion, 
enabling better discrimination between star formation models and clarifying the influence of variable 
accretion on the PLF. 

Subject headings: stars: formation stars: luminosity function, mass function 



1. INTRODUCTION 

Protostars are born in dense, compact molecular cloud 
cores (McKee & Ostriker 2007). During their earliest 
evolution, protostars are both dim and heavily obscured 
by a dusty envelope. It is an unavoidable observational 
reality that the majority of the star formation process 
occurs while protostars are deeply embedded within their 
natal gas. Consequently, high extinction and significant 
radiation reprocessing render the details of protostellar 
evolution, lifetimes, and the accretion process extremely 
uncertain (White et al. 2007). 

Recent surveys of the nearby star- forming regions have 
successfully obtained large statistical samples of young 
protostars with reasonable completeness down to lumi- 
nosities of ^ 0.1 Lq (e.g., Evans et al. 2009; Enoch 
et al. 2009; Dunham et al. 2008). High- resolution mil- 
limeter emission maps tracing dusty envelopes, provide 
a measure of core masses (Enoch et al. 2008). Combined 
with mid to far-infrared data, the available wavelength 
range is sufficient to trace the spectral energy distribu- 
tion (SED), from which the total luminosity may be es- 
timated. 

Using the infrared spectral slope, observed sources can 
be divided into four classes that can be approximately 
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mapped to evolutionary stages (Andre et al. 2000). Class 
protostars are heavily obscured by a dusty envelope, 
such that most of the radiation falls in the sub- millimeter 
band. During the Class I phase the protostar, while 
still embedded, becomes less obscured and may be sur- 
rounded by a thick circumstellar accretion disk. By the 
Class II phase, the now pre-main sequence star has ac- 
creted or expelled most of the initial envelope mass and 
produces little sub-millimeter emission. The remaining 
gas lies in thin accretion disk surrounding the star. Sig- 
natures of outflows may be apparent during both the 
Class I and Class II phases. During Class HI, the disk 
dissipates and the source approaches the main sequence. 

Despite this straightforward picture, cataloging 
sources and definitively mapping them to a physical stage 
remains complicated. Geometric effects shift objects over 
the Class O/I boundary, distorting the correlation be- 
tween physical stage and SED characteristics. Edge-on 
Class II protostars with higher extinction may be mis- 
classified as Class I, while Class I sources viewed down 
the outflow cavity resemble Class II sources (Masunaga & 
Inutsuka 2000; Robitaille et al. 2006). Variability in the 
accretion rate may cause the protostar to oscillate across 
class boundaries (Dunham et al. 2010). Measurements of 
the millimeter emission used as a proxy for the envelope 
mass can be used to distinguish between embedded, i.e. 
Class and Class I objects, and non-embedded. Class II 
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objects (Enoch ct al. 2008), although objects with thick 
disks along the line-of-sight may still be misclassificd. 

During the Class and I phases, luminosity due to 
accretion likely dominates the total radiative output 
(Evans et al. 2009; Dunham et al. 2010). Consequently, 
upper limits for the accretion rates may be inferred from 
the luminosity (Enoch ct al. 2009). This information 
gives clues about the formation timcscalc and the ac- 
cretion process during which the protostars are deeply 
embedded and cannot be directly imaged. However, the 
large distribution of observed luminosities, and in par- 
ticular the significant number of dim protostars, creates 
a picture that is largely inconsistent with predictions of 
some star formation models. 

1.1. The Luminosity Problem 

The accretion rates inferred and the star formation 
times predicted by theoretical models combine to sug- 
gest that protostars are on average too dim (e.g., Kenyon 
et al. 1990; Young & Evans 2005; Enoch et al. 2009). The 
typical accretion luminosity of a protostar 
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exceeds the observed average luminosity of about 2 Lq 
and significantly exceeds the observed median luminosity 
of 0.9 Lq (Enoch et al. 2009). Here we have adopted a 
typical mass accretion rate of m = 2.5 x 10~^A'/q yr^^ 
(McKee & Ostriker 2007) and a typical protostellar ra- 
dius of r = 2.5 Rq (Stabler 1988; Hartmann et al. 1997). 
The factor /acc is the fraction of the gravitational poten- 
tial energy that is radiated away. Note that Lace is the 
total accretion luminosity, including any emission from 
an accretion disk around the star. There are three factors 
that can reduce the accretion efficiency /acc below unity: 
First, some of the energy of the gas that accretes onto 
the protostar can be extracted by a hydromagnetic wind 
(e.g., Ostriker & Shu 1995). If half the energy lost by 
the disk is mechanical instead of radiative, then, since 
half the total potential energy is lost by the disk, this 
would give /acc — | • Second, some of the energy in the 
in-falling gas can be absorbed by the star; Hartmann 
et al. (1997) argue that this is not significant for accre- 
tion rates ^ 10~^Mq yr~^, and this has been confirmed 
by Commergon et al. (2011) so long as the accretion flow 
is transparent to optical radiation. Finally, some of the 
accretion energy is used in dissociating and then ioniz- 
ing the accreting gas (Tan & McKee 2004), but this is 
significant only for very low masses. Thus, of these three 
effects, the mechanical loss of energy is the only one of 
significance. As we shall see in Section 3.4, episodic ac- 
cretion can reduce the observed mean luminosity below 
the true mean and therefore contribute an effective re- 
duction in /acc- 

1.2. Review of Past Work 

A number of mechanisms have been suggested to re- 
solve the luminosity problem. A successful solution must 



account for both the mean luminosity of the distribu- 
tion and the spread of luminosities over several orders of 
magnitude. The first indication of the discrepancy be- 
tween accretion and luminosity was observed in Taurus 
by Kenyon et al. (1990); Kenyon & Hartmann (1995). 
At that time, observations of only a handful of young 
protostars existed, and the authors speculated that the 
problem could be resolved by, among others, short, un- 
observed periods of high accretion or increased ages of T 
Tauri stars. 

Fletcher & Stabler (1994a,b) performed the first 
derivations of protostellar mass and luminosity functions. 
Assuming constant accretion rates and a fixed time pe- 
riod during which stars formed with a constant star for- 
mation rate, they derived the time-dependent luminos- 
ity function for a cluster of evolving protostars. Their 
model included the luminosity contribution from Class I 
to main-sequence stars. For steady-state star formation, 
our constant accretion-rate model produces a protostellar 
protostellar luminosity function similar to that derived 
by (Fletcher & Stabler 1994b) when pre-main sequence 
stars are excluded (sec §3). 

Myers et al. (1998) used a semi-analytic model to de- 
scribe the evolution of the protostar, disk, and envelope 
system. With simple radiative transfer, they generated 
bolometric temperature and luminosity tracks from Class 
to the zero-age main sequence stage. In order to repre- 
sent the decline in accretion luminosities with time, the 
models assumed an exponentially decreasing accretion 
rate. For some parameters, the L — T tracks were able to 
reproduce the mean protostellar luminosities. However, 
the breadth of the distribution and the relative number 
of low-luminosity sources were not addressed. 

Monte-Carlo radiation transfer modeling by Whitney 
et al. (2003) underscored the the dependence of the SED 
on core geometry. In particular, they showed that the 
inferred bolometric luminosity may vary by 50% from 
the true luminosity depending upon the orientation of 
the disk and outfiow cavity relative; to the line-of-sight. 
For a sufficient number of observed protostars sampling 
all geometric projections, the orientation can only alter 
the distribution of the luminosities, not the mean. 

Young & Evans (2005) performed 1-D radiative trans- 
fer observations of the inside-out collapse model devel- 
oped by Shu (1977). These calculations improved upon 
the Myers et al. (1998) work by following the hydrody- 
namic evolution of the protostellar core. They included 
six contributions to the total luminosity with the result 
that their L — T tracks agreed only with the brightest 
protostars, further illuminating the discrepancy between 
observed luminosities and predictions of low-mass star 
formation. 

Recent extensive surveys of five local star-forming re- 
gions have provided more comprehensive statistics and 
observational data of the earliest stages of protostellar 
evolution (Enoch et al. 2009; Evans et al. 2009) This 
work demonstrated that Class and Class I sources have 
comparable mean luminosities, suggesting that proto- 
stars may accrete significantly for longer periods of time 
than previously assumed. Extending the accretion epoch 
from 0.1 Myr, as first assumed by Kenyon et al. (1990), to 
0.5 Myr presents one possible solution to the luminosity 
problem. 

Dunham et al. (2010) revisited the Young & Evans 
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(2005) methods and explored several improvements, in- 
cluding updated opacities, 2-D effects such as disks, core 
rotation, outflow cavities, and variable accretion. With 
all the additions, although predominantly due to the in- 
clusion of episodic accretion, good statistical agreement 
between the observed and modeled temperature and lu- 
minosity distributions was possible. They modeled vari- 
ability by halting accretion from the disk onto the pro- 
tostar until the disk mass reached 20% of the star mass, 
whereupon the star was allowed to accrete at a constant 
rate of 10~* Mq yr~^ until the disk mass was exhausted. 
Since other terms contributing to the total luminosity 
(e.g., accretion from the envelope onto the disk and the 
stellar luminosity) are small compared to the accretion 
luminosity, the model assumes that protostars exist in a 
very low luminosity state for most of their formation time 
and accrete most of their mass during bursts. Statisti- 
cally, their models predict that protostars must spend 
~ 1% of the time radiating at L > IOOLq. This is 
only marginally consistent with the observed sample of 
112 protostars, which contains no sources more luminous 
than 76 Lq . The models also have a star formation time 
< 0.2Myr, which is less than the star formation time 
inferred by Evans et al. (2009). Future surveys contain- 
ing more objects are necessary to test the importance of 
episodic accretion. 

There is debate whether allowing for episodic accretion 
in evolutionary models can reproduce the dispersion of 
low-luminosity sources on the Hertzsprung-Russell dia- 
gram and thus account for a portion of the age spread 
inferred for members in young, low-mass star clusters 
(Baraffe et al. 2009; Hosokawaet al. 2011). If so, episodic 
accretion could also explain lithium depletion measure- 
ments in young stars (Baraffe & Chabrier 2010). How- 
ever, this solution is valid only if nearly all the accretion 
energy is efficiently radiated away from a small fraction 
of the stellar surface. Other mechanisms such as varying 
initial conditions and assumptions about the radiative 
properties of the accretion flow, which are not well con- 
strained, may also contribute to an apparent age spread 
(Hosokawa et al. 2011). 

1.3. Model Assumptions 

In this paper we derive the Protostellar Luminosity 
Function (PLF) predicted by various star formation mod- 
els and compare with observations (Enoch et al. 2009). 
Comparison with the mean observed luminosity deter- 
mines the mean star formation time; comparison with 
the shape of the PLF (ratio of median to mean and the 
standard deviation) tests the validity of the models. Our 
basic formalism (McKee & Offner 2010a, henceforth Pa- 
per I) can be adapted to consider a time- varying rate of 
star formation, bursty accretion, and an arbitrary stellar 
initial mass function. However, in this paper we base our 
comparison on four main assumptions. First, we treat 
only the cases in which the star formation rate is either 
constant or smoothly accelerating in time (e.g. Palla & 
Stahler 2000). Second, we assume the accretion rate is 
a smooth function of time during the accretion phase; in 
particular, we assume it can be expressed as a function 
of the current mass and the final mass of the protostar. 
Rather than excluding stochastic variability completely, 
we assume that any time-varying component is statis- 
tically rare in the samples with which we compare and 



therefore unlikely to be included the data (sec §5.3 for 
a detailed discussion). Third, we assume that the accre- 
tion rate onto the protostar, which can be inferred from 
the protostellar luminosity, tracks the accretion rate onto 
the protostar-disk system over the majority of the life- 
time of the protostar. Finally, we adopt an individual- 
star Chabrier (2005) IMF truncated at an upper mass 
limit m„. 

Through our derived PLF models, we investigate the 
luminosity problem in the context of different star forma- 
tion theories and explore variations such as an accelerat- 
ing star formation rate and accretion rates that taper off 
over time. In §2 we review the Protostellar Mass Fimc- 
tion, which is described in detail in Paper I. In §3 we 
derive the PLF and relative statistics. We compare with 
observations of local star-forming regions in §4 and dis- 
cuss our results in §5. Approximate analytic results for 
the PLF for different star formation models are given in 
the Appendix. 

2. THE PROTOSTELLAR MASS FUNCTION (PMF) 

In Paper I we determined the mass function of pro- 
tostars (the PMF) in terms of the IMF and the accre- 
tion history of the protostars. We assumed that the ac- 
cretion history could be expressed in terms of the cur- 
rent protostellar mass, m, and its final, stellar mass, 
TO/. Henceforth, to and toj are measured in units of 
solar masses and to is in units of solar masses per year. 
We expressed the IMF as tp{mf), where ip{mf)dlnmf is 
the fraction of stars born with final masses in the range 
druf. We then defined the bi-variate PMF, 7/;p2, such 
that 'ipp2 d In TO d In to,/ is the fraction of protostars in a 
region of star formation with masses in the range dm 
and final masses in the range dm/. We showed that for 
steady (i.e., non-accelerating) star formation 

Mm,mf) = —^, (3) 

where t f is the formation time for a star of mass to / and 
{t f) is the average of t/ over the IMF. The PMF, ipp{m), 
is just the integral of the bi-variate PMF, 

tpp{m) = / tpp2{m,mf)dlnmf, (4) 

where the lower limit of integration is 

mf^i = max(TOf, to) (5) 

and the most and least massive stars formed by the clus- 
ter are to„ and mf , respectively. 

To determine the PMF, we required the accretion his- 
tories of the protostars. As noted above, we assume that 
these accretion histories, apply to the gas reaching the 
protostellar surface as well as to the gas that accretes 
from the ambient medium onto the protostar-disk sys- 
tem. We considered several different models: 

i) Inside-out collapse of an isothermal sphere (IS, Shu 
1977), 

rh = rhis = 1.54 x 10-^(T/10 Kf^'^. (6) 

h) The Turbulent Core model (TC, McKee & Tan 
2002, 2003), in which the initial core is presumed 
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to be supported by turbulent niotions instead of 
thermal pressure, 



where 



m = 3.6 X lO-^E^Y" (^-^y m/^^ 
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where m is in Mq yr~^, Ed is the surface density of 
the clump in g cm^^ in which the stars arc forming, 
and the parameter j is determined by the density 
profile of the core; we follow McKee & Tan (2003) 
in setting i = 5. Tan & McKee (2004) suggested 
that the coefficient in the accretion rate in equa- 
tion (8) be increased by a factor ^ 2.6 in order to 
allow for infall at the beginning of the inside-out 
collapse. We retain the value given by McKee & 
Tan (2003), but note that the overall normaliza- 
tion of the accretion rate remains uncertain. 

iii) A simplified model to represent competitive accre- 
tion (CA, Bonnell et al. 1997, 2001), 



2/3 

m = mi\ — ] m/ 



(9) 



where mi is the accretion rate for m = m/ = 1. 
Note that this accretion rate has the property that 
stars of all masses form in the same amount of time. 

All these forms for the accretion history can be sum- 
marized in the expression 



m 



\mfj 



where the model exponents are given by: 
IS: j=j>=0 

TC:j = -,j/ = - 
J 2'-" 4 



1. 



(10) 

(11) 
(12) 

(13) 



We also considered a two-component turbulent core 
model (2CTC), a blend of the IS and TC models: 

_,2 m \ 
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where TZm = mxc/'Ttis; wc adopt TZm = 3.6 as in Pa- 
per I. This model is similar to the TNT model of Myers 
& Puller (1992). There is some evidence that a two- 
component accretion model may be more physical in the 
CA case as weU (i.e., a 2CCA model). Smith et al. (2009) 
show that, at least initially, competitively accreting pro- 
tostars are surrounded by a small bound envelope. This 
suggests that at early times the protostars may accrete 
via a Shu-like constant accretion rate. However, the au- 
thors find that the envelope mass is not well correlated 
with the final mass of the star, which they suggest indi- 
cates that the CA phase dominates. For comparison, we 
define a 2CCA model with: 
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According to Smith et al. (2009), accretion of the core en- 
velope should contribute less than half of the final mass. 
Therefore, we adopt T^rri, CA = 3.2, which is determined 
by assuming that a star of average mass accretes half its 
mass from its envelope. 

In the fiducial cases above, the accretion rates increase 
monotonically over the protostcllar lifetime. For the 
single-component models, the time to form a star, tf is 
a power-law function of the final mass. 



tf = tfiUlf 



where 



1 



(17) 



(18) 



(1 - j)mi 

is the time to form a star with m = 1. In reality, we 
expect accretion rates to gradually taper off (e.g., Myers 
et al. 1998; Myers 2010), at least for the IS and TC mod- 
els. Tapered accretion is supported observationally since 
Class I sources are not significantly brighter than Class 
sources, and in some cases, the most luminous young 
source is actually Class (e.g., Evans et al. 2009). Class I 
sources are also twice as numerous, suggesting that high 
accretion rates, which depend upon the presence of dense 
infalling gas, cannot be sustained throughout the major- 
ity of the formation time (Enoch et al. 2009). In com- 
petitive accretion models, the decline in accretion rates is 
assumed to be abrupt, so untapered models are likely to 
be the best approximation for the CA case. If we assume 
that the decline is linear in time, the tapered accretion 
case can be written: 

m = mo(m, m/) 



1 



(19) 



where mo (m, mf) is the untapered accretion rate for m = 
mf — 1 Mq. We can then express the accretion rate for 
both the tapered and untapered cases as 

nl 1/2 



m = mo 1 m r 

' * m/ ' 
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(20) 
mf = 

(21) 



where mo, 1 is the untapered accretion rate for m ■ 
1 and 

untapered, n = 0, 
tapered, n = 1, 

The star formation time for single-component models can 
then written as: 

tf = tfimf^-^f{l + Sni). (22) 
In the case of accelerating star formation we assume a 
birthrate that increases exponentially with time, i.e., 

K oc e(*-*'")/^ (23) 

where tm is the age of a star with mass m and final mass 
mf and t is the characteristic acceleration time. For 
an accelerating star formation rate, the mean formation 
time is no longer equal to the observed star formation 
time. Instead, 

(t„)(l-e-*//-) 

(*/)obs ^ 



(e*//-^)(l - e-'"/^) 
{tu) 



t(1 - e-^*")/-^) 



(24) 
(25) 
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where tu is the Hfetime of Class II sources. For tu = 2 
Myr and r = 1 Myr, the actual mean star formation time 
of the accelerating cases, {tf), is a factor ~ 2.3 smaller 
than the value inferred from observations of the numbers 
of protostars and young stellar objects, (i/)^^^- 

3. THE PROTOSTELLAR LUMINOSITY FUNCTION (PLF) 
3.1. Definition 

Our objective is to determine the protostellar luminos- 
ity function, ^'p(L), where 'i>p{L)d\n L is the fraction of 
protostars that have luminosities in the range dL. The 
protostellar mass and accretion rate are related to the 
accretion luminosity through equation (1), which gives 
L oc mm/r. Unfortunately this relation is complicated 
by the fact that the protostellar radius is a function of 
both the mass and the accretion rate, r = r(rn, rh), which 
does not have a simple analytic representation (Stabler 
1988; Hartmann et al. 1997). We use the routine de- 
scribed by OfFner et al. (2009), which agrees with the 
results of Stabler (1988), to determine the protostellar ra- 
dius (see §3.3). This model exhibits evolution similar to 
the model with shock boundary conditions in Hosokawa 
ct al. (2011). Consequently, the subsequent evolution of 
the stars will have ages that are consistent with previ- 
ously calculated non-accreting isochrones. 

We proceed by using equation (10) to determine 
r(m,m/) from r(m,m), and then solving the accretion 
luminosity equation (1) for to(L, my ). Since there are 
two independent variables, we use bi-variate distribution 
functions: The fraction of protostars in the luminosity 
range dL and final mass range druf is the same as the 
fraction protostars with masses in the range dm and final 
mass range drnf. 

'i!p2{L,mf)dlnLdlnmf = il)p2{m,mf)dlnm dlnnif , 

(26) 

where 'ipp2 was determined in Paper I (see eq. 3). The 
PLF is then obtained by integrating '^p2{L,mf) over 
m/. 
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In the Isothermal Sphere case, the accretion rates are 
independent of the final mass. Consequently, the PLF 
reduces to: 

ML) r".,„^ V'(m/) 



din TO/ 



1 



In r I 
d\nm\ 



(29) 



3.2. The Mean and Standard Deviation of the 
Luminosity 

In general, the mean luminosity can be evaluated as 
(L) = J L'^p{L)d\nL. However, it is more instructive to 
return to the bi-variate luminosity function and define 
the mean accretion luminosity: (26): 
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where equation (26) enabled us to eliminate the depen- 
dence on TO for steady star formation and where 



1 



r(TO/) fnf^ Jq r{m,mf) 



mdm 



(32) 



is the mass-averaged harmonic mean radius. Note that 
equation (31) is only valid in the case of steady star for- 
mation, where the accretion luminosity is proportional to 
the binding energy of the stars formed. In order to eval- 
uate the mean luminosity, we take advantage of the fact 
that the stellar radius depends primarily on the instan- 
taneous mass TO and only weakly on the final mass to/ 
in order to define a mean radius in terms of the current 



mass m. 



r{m) 



mfdvif 
r{m,mf) 



(33) 



This approximation is reasonably accurate for the value 

of the upper mass limit we consider here (m„ = 3); for 
larger values of m„, it would be preferable to include a 
weighting with the IMF. 

For protostars with to/ ^ 1 Mq, the stellar luminosity, 
i.e., the luminosity due to nuclear burning and Kelvin- 
Helmholtz contraction, makes only a small contribution 
to the total luminosity, so that (L) ~ (iacc)- However, 
for more evolved protostars beginning late in the Class I 
stage, accretion diminishes and no longer dominates the 
total luminosity (White & Hillenbrand 2004). For the 
young, low-mass, embedded sources we consider here, 
the accretion luminosity is much greater than the nu- 
clear luminosity provided m I.^Mq. For such stars 
with tapered accretion, the accretion luminosity dom- 
inates during the majority of the evolution, i.e., while 
TO < 0.95 To/. For stars with to > 1.5Mq, the nuclear 
luminosity becomes important and the accretion lumi- 
nosity is not a good approximation of the total. However, 
these stars comprise only a small part of the total sam- 
ple. For example, the maximum expected protostellar 
mass, TOinax = 2.3Mq, for a cluster with to„ = 3Mq and 
J\fp = 112 (see Figure 7 in Paper I) in the CA case corre- 
sponds to a maximum ratio of Lacc/^tot = 0.65. In the 
fiducial IS case, Wmax = 2.9M0 yields iacc/-^tot = 0.16, 
which reflects both the lower accretion rate of the IS 
model and the higher maximum mass for a cluster with 
a fixed number of protostars in this model. 

The standard deviation is a useful metric for character- 
izing the breadth of the protostellar luminosity distribu- 
tion. Since the luminosity spread may encompass several 
ord(Ts of magnitude, we adopt the standard deviation of 
the logarithm of the protostellar luminosity: 



o-(logL) : 



j (logL)2*p(L)dlnL 
J \ogL *p(L)dlnL 



1 1/2 



(34) 



This metric has the advantage that it is dimensionless 
and, like the mean, it can be computed solely in terms of 
the protostellar mass function if L{m,mf) is provided. 
The ratio of the median to the mean, also a useful di- 
mensionless number, measures the skewness of the dis- 
tribution. 
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Fig. 1. — Stellar radius versus mass for a constant accretion rate 
of IQ-^MQyr-i (solid), the IS model (m = 9 X 10-^M©yr-l, 
dotted), the TC model (dashed), and the CA model (dot-dashed) 
where {tf) = 0.44 Myr and nif = 10 Mq. For comparison, mod- 
els calculated by Hosokawa & Omukai (2009) for constant ac- 
cretion rates of 10~^ Mq yr~^ and 10~^MQyr~^ are shown by 
the plus signs and diamonds, respectively. Stars indicate the 
Palla & Stahler (1991) model for a constant accretion rate of 
IQ-^Mgyr-i. 

3.3. Stellar Evolution Model 

In order to determine the PLF predicted by a given 
theory of star formation, we must first construct a com- 
prehensive luminosity model that takes into account both 
stellar evolution and the accretion history. We adopt the 
one-zone protostellar evolution model developed by Tan 
& McKee (2004), which is described in detail in OfFner 
et al. (2009). This model has been calibrated to repli- 
cate the more detailed calculations of Palla & Stahler 
(1991, 1992) and updated to reflect the recent recom- 
mendations of Hosokawa & Omukai (2009). It corre- 
sponds to the "hot" stellar evolution model described 
in Hosokawa et al. (2011), where an accretion flow di- 
rectly hits the stellar surface and forms an accretion 
shock front. Although the gas may be channeled through 
a disk, the model assumes that it is sufHciently thick 
such that the accretion column covers most of the stellar 
surface. In practice, the initial physical conditions and 
accretion flow properties are not well known, introduc- 
ing uncertainties in the stellar radius of a factor of ~ 3 
for m < 0.2Mq. (e.g., Hosokawa et al. 2011). Figure 
1 shows the stellar radius predicted by this model as a 
function of stellar mass for different accretion histories. 
For m = lO~^M0yr~^, the model has a maximum dis- 
agreement with Hosokawa & Omukai (2009) of ~ 20% 
when m < 3Mq and generally agrees within 5%. The 
model has a maximum disagreement of a factor of ^2 
at m = 1.6Mq for m = lO~^M0yr~^, but it gener- 
ally is within a factor of 1.2 Hosokawa & Omukai (2009). 
Although our sub-grid model has been thoroughly de- 
scribed in Offner et al. (2009), we give a brief summary 
here. 

The model treats the protostar as a polytrope and fol- 
lows the stellar contraction by enforcing energy conser- 
vation. The model is characterized by six stages, cul- 
minating in the arrival of the protostar on the zero age 
main sequence. In the "pre-collapse" stage, the collapsing 



gas densities are sub-stellar. In the second "no burning 
stage", the densities are sufficient for the dissociation 
of H2, but are too low for deuterium burning. In the 
following "core deuterium stage at flxed Tc," deuterium 
burning ignites in the core. Once the supply of deu- 
terium is depleted, the core temperature rises and the 
protostar enters the "core deuterium burning at variable 
Tc." Eventually high temperatures in the core termi- 
nate convection and halt the flow of deuterium into the 
core, whereafter the protostar enters the "shell deuterium 
burning" stage. Finally, the central temperature reaches 
10'' K, and the star moves onto the main sequence. In 
calculating the total luminous output, the zero- age main 
sequence luminosity serves as an approximation of the 
interior luminosity, which results from the combination 
of deuterium burning and Kelvin-Helmholz contraction. 
We evolve this model in combination with the accretion 
rates specifled by the star formation models and generate 
luminosity tables as a function of the instantaneous and 
final stellar masses for each case with both untapered 
and tapered accretion. For an accelerating star forma- 
tion rate, we assume that the physical parameters are set 
by the initial conditions of the collapse and hence do not 
themselves accelerate. Thus, for the untapered acceler- 
ating and tapered accelerating star formation cases, we 
use the untapered and tapered luminosity tables, respec- 
tively. 

We use these tabulated values of L(m, m/) to obtain 
the mean and standard deviation of the luminosity di- 
rectly. However, in the model, L and r undergo discon- 
tinuous jumps as the deuterium state changes. This is 
problematic in calculating the PLF, which requires the 
derivative of L. To circumvent this issue, we use a poly- 
nomial fit to equation (33). We re- normalize the result 
using the directly derived value of (L). This strategy 
simplifies the integral while preserving the trends in the 
sub-grid model and eliminating the weak dependence of 
L on mf. 

3.4. Episodic Accretion 

Observations of high-luminosity variable sources, such 
as the prototype FU Ori, suggest that protostars undergo 
periods of high accretion (e.g., Hartmann & Kenyon 
1996). In the extreme case, protostars may spend most 
of their life in a low-luminosity, low-accretion phase and 
accrete nearly all their mass during short, intense ac- 
cretion bursts. However, only a total of 18 bursting 
sources have been identified within 1 kpc of the Sun 
(Greene et al. 2008). These sources have luminosities 
in the range 20-550 Lq, corresponding to accretion rates 
of 10"^ - 10""* Mq yr-^ The total mass that can be 
accreted in such bursts is limited by the known star for- 
mation rate. Following McKee & Offner (2010a), let /epi 
be the fraction of mass accreted during episodic accretion 
periods: 

f _ ™cpi^^cpi /oe,\ 

where Aicpi is the total time spent in the bursting state 
and (nif) = 0.5 is the mean stellar mass for a typical 
IMF. The time spent in the high accretion state can be 
expressed as: 

Atcpi - (36) 
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where A/'p.cpi — 20 is the number of protostars accreting 
in this state within 1 kpc of the Sun (Greene et al. 2008) 
and the star formation rate, AC, is 0.016 stars/yr within 
1 kpc (Fuchs et aL 2009). This gives Aicpi ^ 1200 yr and 
/cpi ~ 0.25 for an accretion rate of 10~^ Mq yr~^. If the 
protostehar hfetime is 0.5 Myr, then the duty cycle must 
be - 0.2%. 

This estimate suggests that our sample of 112 proto- 
stars is unlikely to contain any bursting sources. For our 
observational comparison, we adopt an episodic accre- 
tion factor, /cpi — 0.25, to reflect the amount of mass 
accreted during outbursts. If no protostars are currently 
undergoing bursts, then the observed sample will have a 
mean luminosity that is 75% of the true time-averaged 
mean. To account for this, we reduce our model mean 
luminosity by a factor of 1 — /cpi = 0.75. In other words, 
we adopt an effective value for /acc in equation (1) of 

/acc,cff = /acc(l " /cpi) = 0.56. (37) 

Note that the Dunham et al. (2010) episodic accretion 



model corresponds to /< 



epi 



0.7 for a star with 



0.35. For /acc = 0.75, this corresponds to /acc, off = 0.23. 
An alternative way of correcting for episodic accretion 
assumes that the protostars accrete at their normal rate 
for a time t non-episodic and then accrete via bursts for a 
brief additional time. In this case, the luminosity of the 
protostars would be unchanged and instead the forma- 
tion time would be correspondingly shorter. We adopt 
the former definition of /epi since it seems more plausi- 
ble and is consistent with simulations exhibiting episodic 
accretion (e.g., Vorobyov & Basu 2005). 

4. COMPARISON WITH OBSERVATIONS 

4.1. Overview of the Observational Data 

Throughout this section, we compare to the Class 
and Class I source luminosities of Serpens, Ophiuchus, 
Perseus, Lupus and Chameleon from Enoch et al. (2009) 
and Evans et al. (2009). This is the same sample used 
in Dunham et al. (2010). It has been corrected for lo- 
calized extinction and carefully analyzed to remove non- 
protostars and non-embedded sources (Mike Dunham, 
private communication). The data are comprised of 39 
Class and 73 Class I objects for a total sample size of 
112. This includes only two sources from Chameleon II 
and one source from Lupus, such that the sample essen- 
tially reflects the properties of Serpens, Ophiuchus, and 
Perseus. The Class and Class I sources have a 0.50 like- 
lihood of being drawn from the same population, which 
suggests that the two groups are quite similar. It is un- 
clear if this is due to similar accretion rates or source 
mixing resulting from inclination effects (e.g., Robitaille 
et al. 2006). 

As discussed by Enoch et al. (2009) and Dunham et al. 
(2008), the bolometric luminosities have an uncertainty 
of ^ 50% due to saturation of the 160 /im band. Enoch 
et al. (2009) estimate that ^ 10% error arises from the 
SED fitting process and 15-25% due to finite sampling 
errors. Although the data we use have improved 350 
jim observations and more careful source selection, un- 
certainty arises from a number of factors and is not well 
constrained. Following Enoch et al. (2009), we adopt 
a 50% upper uncertainty in the bolometric luminosities, 
and we use the non-extinction corrected bolometric lu- 
minosities, which are naturally underestimates, to set 




10.00 



Fig. 2. — The PLF, ^p{L), of the observed protostellar lumi- 
nosities for Perseus, Serpens, and Ophiuchus (Evans et al. 2009). 

a lower error limit (Melissa Enoch and Mike Dunham, 
private communication) . The extinction-corrected mean 
and median bolometric luminosities are then: 



(-^Obs) — 5.3^^'g Lq, 



_i c:+0.7 
^obs,mcd J--^ 



0.4 ^Q- 



(38) 
(39) 

Note that the mean is very sensitive to the most lumi- 
nous sources. The extinction corrected sample of Enoch 
et al. (2009), which is derived from a nearly identical 
source list, has a mean of 3.5 Lq, which falls within the 
error of our adopted mean. The disagreement arises 
from a combination of improved 350 /im data used in 
the Evans et al. (2009) analysis and different treatments 
of the IRAS fluxes. In either case, the luminosity of a few 
of the brightest sources may have been overestimated by 
as much as a factor of two if the 24/xm flux was omitted 
due to saturation (Mike Dunham, private communica- 
tion). Consequently, the upper uncertainty should be 
considered an upper limit on the luminosity rather than 
a one-sigma error estimate. Since the standard deviation 
is also sensitive to the distribution outliers, we use the 
log of the standard deviation for comparison: 



a(logLobs) = Q.lZli dex 



(40) 

A second useful non-dimensional quantity is the ratio of 
the median luminosity to the mean: 



imcd/(i) =0.3 



-1-0.2 
0.1- 



(41) 



However, the associated uncertainties of this ratio are 
somewhat uncertain since the median and mean errors 
are correlated. 

In this sample, the mean luminosities of the individ- 
ual observed star forming regions and their shapes are 
similar, even though the number of sources in each is 
different. The PLFs of Perseus, Serpens, and Ophiuchus 
are shown in Figure 2, where the regions have mean lu- 
minosities 5.2 Lq, 5.0 Lq, and 6.2 Lq, respectively. A 
Kolmogorov-Smirnov (K-S) test indicates that Perseus 
and Serpens have a 0.39 likelihood of of being drawn 
from the same parent distribution, while Perseus and 
Ophiuchus have a 0.45 likelihood and Ophiuchus and 
Serpens have a 0.86 likelihood. This high level of agree- 
ment, despite the apparent dissimilarity of the distribu- 
tions shown by Figure 2, is partially due to the small 



8 



Offner & McKee 



number of protostars 20) in Serpens and Ophiuchus. 
While comparing our models with the individual regions 
would be ideal, small statistics requires that we use the 
total protostellar distribution over all the clouds. 

To derive the model distribution, we require the mass 
of the most massive star forming in the cluster, m„. The 
brightest Class II source in Evans et al. (2009) suggests 
an upper mass of ~ 3 Mq. Since the number of Class II 
sources is several times larger than the sample of Class 
I and Class sources, this should statistically be a good 
upper limit. However, statistical variations or changes in 
the star formation rate may impact the upper mass of 
the currently forming stars. Inspection of the most lumi- 
nous Class I source also suggests an upper mass of 3Mq, 
assuming that its luminosity is dominated by stellar ra- 
diation rather than accretion. This latter estimate is also 
very uncertain since we can't discount that the source is 
undergoing an outburst and its luminosity is dominated 
by accretion. However, we can be fairly confident that 
m„ = 3 Mq is a good upper limit. 

The accretion timescale of the models is constrained 
by the observed protostellar lifetime, i.e., how long pro- 
tostars spend in the main accretion phase. Both Stage 
(Menv > m) and Stage I (m > Mgnv > O.IM©) pro- 
tostars experience significant accretion, since a signifi- 
cant fraction of the total gas is contained in the envelope 
(Crapsi et al. 2008). These stages roughly correspond to 
Class and Class I provided that the Class I envelope 
mass exceeds 0.1 Mq. Enoch et al. (2009) estimated 0.1 
±0.02 Myr for the Class lifetime. This is longer than 
the Class lifetime of 1 — 3 x 10^ yr reported by Andre 
et al. (2000), because these authors base their lifetime 
solely on data from Ophiuchus, which is now recognized 
to have a deficit of Class objects compared to other 
star forming regions (Enoch et al. (2009), who also in- 
clude one more Class object than Andre et al. (2000) in 
their Ophiuchus sample.) For local star forming regions, 
Evans et al. (2009) report an average combined Class 
and I lifetime of 0.44 Myr. The uncertainty in this esti- 
mate arises in part from statistics and Class I/II source 
confusion, but it is dominated by the uncertainty in the 
disk lifetime (2 ± 1 Myr), which is necessary to calibrate 
the ages. Altogether, the mean protostellar lifetime is 
(tf) = 0.44 ± 0.22 Myr in the Evans et al. (2009) sam- 
ple. Their estimates of the mean protostellar lifetimes in 
individual clouds, Ophiuchus (0.30 Myr), Serpens (0.46 
Myr), and Perseus (0.62 Myr), are within the errors of 
the overall mean. This lifetime is significantly longer 
than previous estimates, which have adopted cither the 
Class lifetime (e.g., Kenyon et al. 1990) or the core free- 
fall time (~0.1 Myr for a 0.5Mq star; e.g., Myers et al. 
1998; Young & Evans 2005). Shorter lifetimes exacer- 
bate the luminosity problem. As discussed by McKce & 
Offner (2010a), one possible solution is "slow accretion," 
in which accretion rates are reduced by protostellar out- 
flows or lengthened disk lifetimes. Estimated protostellar 
lifetimes on the order of 0.5 Myr lend credence to a slow 
accretion picture. 

Some studies of Class I protostars have found that 
stellar luminosity dominates the total luminosity, sug- 
gesting that accretion has already diminished Muzerolle 
et al. (1998); White & Hillenbrand (2004); Connelley & 
Greene (2010). However, in-depth study often reveals 
that sources classifed as Class I are actually edge-on or 



heavily obscured Class II sources (White & Hillenbrand 
2004; van Kempen et al. 2009; Heiderman et al. 2010). 
The observational sample here uses a minimum gas mass 
criteria of O.IM© for inclusion in the sample, which is an 
indicator that accretion is still underway. However, the 
authors do not probe for tracers of dense, warm gas, such 
as emission from higher transtions of HCO+ and C'^^O, 
which more conclusively separates embedded protostars 
from obscured Class lis van Kempen et al. (2009); Hei- 
derman et al. (2010). 

4.2. The Mean Luminosity, (L) 

The mean of the luminosity distribution serves as a 
simple one-dimensional statistic with which to compare 
the models and observations directly. In Figure 3, we plot 
the mean luminosity as a function of the most massive 
star forming in the cluster for each of the accretion mod- 
els. We truncate the IS model at 5Mq since it is unreal- 
istic for describing high-mass star formation. In Table 1, 
we give the model values calculated assuming complete- 
ness of the extinction-corrected data down to 0.05^© for 
m„ = 3Mq. We plot these tabulated values and the 
observed mean ((Lobs) = 5.3l^;9 Lq) from Evans et al. 
(2009) in an inset plot in Figure 3. As with the mean 
protostellar mass derived in Paper I, the mean luminos- 
ity increases strongly as a function of the upper stellar 
mass. Unfortunately, for m„ = 3Mq, the model spread 
is small, making it difficult to discriminate between mod- 
els. In contrast, the model means diverge significantly for 
clusters with large m„ in all cases, suggesting that more 
complete observations of high-mass star-forming regions 
would be useful to distinguish models based solely upon 
the mean luminosity. Note that the model means depend 
upon the uncertain star-formation timescale, so that it is 
not possible to make an independent comparison of the 
models and observations (see discussion in §5). However, 
better constraints on the star formation time in the fu- 
ture should increase the discriminating value of the mean 
luminosity in comparing models. 

As shown in the plot insets, the mean observed lumi- 
nosity falls above the models in all the non-accelerating 
cases cases. The untapered TC and CA means, IS ta- 
pered mean, and all accelerating means are consistent 
with the observational error. ^ This clearly indicates that 
there is no luminosity problem in the traditional sense. 
The resolution is a result of the longer protostellar life- 
time adopted from Evans et al. (2009) and an effective 
accretion efficiency, /acc,eflf = 0.56 due to a radiative ef- 
flciency of 75% and allowance for episodic accretion at 
the level of 25%. Altogether this reduces the predicted 
luminosities for the non-accelerating cases by a factor of 
- 3. 

The mean luminosities in the accelerating cases are 
up to 30% lower than in the fiducial non-accelerating 
cases for a fixed value of {tf) because the accelerating 
cases have more low-mass protostars. However, for a 

given observed value of {tj}^^^, the mean luminosities 
are raised since their formation time is a factor ~ 2.3 

^ The CA luminosity distribution extends to both the highest 
and the lowest luminosities, so that it yields the largest mean when 
a cutoff of O.O5-L0 is applied and the curve is re-normalized thus 
weighting the highest luminosities more heavily; the IS case is least 
affected by this luminosity truncation since it is strongly peaked 
around a luminosity much higher than the cutoff. 
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TABLE 1 

Mean luminosity {L){Lq) for {tf)^^^^ = 0.44 Myr; (Lobs) = 5.3 



+2.6 
-1.9 





Non- Accelerating 


Accelerating " 


Model 


Untapered 


Tapered*" 


Untapered 


Tapered*" 


Isothermal Sphere (IS) 


3.27 


3.81 


4.50 


4.49 


2C Turbulent Core (2CTC)'^ 


2.97 


3.17 


5.25 


5.44 


Turbulent Core (TC) 


3.63 


3.09 


5.70 


4.44 


2C Competitive Accretion (2CCA)<* 


2.91 


2.68 


4.97 


4.14 


Competitive Accretion (CA) 


4.26 


3.26 


7.05 


4.73 



"n = 



1 Myr; (tu) 
1 

= 3.6 
, CA = 3.2 



: 2 Myr; {tf) ~ {tf),^j2.3 = 0.24 Myr. 



smaller, as shown in equation (25). It is this adjustment 
that accounts for the good agreement of the accelerating 
models with the observations. 

Allowing the accretion rates to taper off during the 
accretion period has a varying effect on the mean lumi- 
nosities. In the IS case, the accretion rate is no longer 
independent of mass, increasing the mean luminosity. In 
comparison, the other models accrete at higher rates at 
early times, lowering the mean luminosity. Adding taper- 
ing also differentiates the models slightly at m„. How- 
ever, when the curves are re-normalized using the ob- 
servational completeness limit, the differences are mini- 
mized and in the TC and CA cases the mean increases. 
Among the tapered, non-accelerating accretion models, 
only the IS model falls within error under the constraint 
that the mean star formation time is 0.44 Myr. 

4.3. Median Luminosity 

While the mean is sensitive to the maximum luminos- 
ity, which observationally may be prone to both over- 
estimation and statistical fluctuations, the median serves 
as a better proxy for the peak of the distribution. In 
this case, systematic errors in the spectra fitting and ex- 
tinction corrections are likely to dominate the error in 
the observed value. Table 2 gives the medians for each 
of the cases for m„ = 3Mq and {t/)^^^ = 0.44 Myr, 



where the observed median is 1. 



For the fidu- 



cial case, only the 2CTC and TC models are within error. 
When tapered accretion is included, all but the IS model 
agree. For an acclerating star formation rate, the unta- 
pered 2CTC, 2CTC, TC and tapered CA and TC models 
are within error. 

4.4. The Star Formation Timescale, {tf) 

We see that in some cases the mean and median lumi- 
nosities are quite different from the observed values when 
we fix the star-formation time at {tf)^^^ = 0.44 Myr. 
However, this time scale is uncertain by a factor 2. 
Henceforth, rather than fixing the observed lifetime, we 
derive the average formation time, (tf), by adjusting the 
models such that in all cases the average model luminos- 
ity agrees with the average observed luminosity, (L) = 
(Lohs)- Figure 4 shows the star formation timescale 
derived from the observed mean luminosity versus the 
mean model luminosity obtained above from the obser- 
vational timescale reported by Evans et al. (2009). The 
timescales are also listed in Table 3, which includes the 
dimensionless parameter values for each model. For non- 
accelerating models, the mean formation time is the same 



as the value that would be inferred from observation by 
comparing the number of protostars with the number of 
Class II sources. For accelerating models, however, the 
formation time is significantly less, {tf) ~ (t^)^j^_,/2.3, for 
the parameters we have adopted; both values are listed 
in the table. For models with non-accelerating accretion, 
the figure strongly suggests that the actual formation 
time, during which the majority of a protostar's mass is 
accreted, is 



(t/)=0.3±0.1 Myr. 



(42) 



This is longer than the observed Class lifetime of 0.1 
Myr (e.g., Enoch et al. 2009), which implies that signifi- 
cant accretion continues into the Class I phase. 

Alternatively, one could use the median rather than 
the mean to adjust the models. This would result in 
an inferred distribution of {tf) more evenly distributed 
above and below {tf)^^^. For example, the untapered, 
nonaccclcrating CA case would require a time half the 
length of the observed formation time, while the unta- 
pered, nonaccelerating IS case would require a time twice 
as large. However, adjusting the mean gives better agree- 
ment between the overall distributions (see Section 4.5). 

The vertical error bars on {tf)^^^^ arise from the uncer- 
tainty in the adopted disk lifetime: 2±1 Myr. Another, 
smaller source of error arises from the possible misclassi- 
fication of sources. In comparing the observed formation 
time with that predicted by the models summarized in 
Table 3, further uncertainty is introduced by the uncer- 
tainty in the mean luminosity of the sample, which has 
been used to normalize the models. Models without ac- 
celeration have an average formation time of about 0.3 
Myr, which would require disk lifetimes close to 1 Myr, 
smaller than most of the results cited by Evans et al. 
(2009). Models with an acceleration time r = 1 Myr 
have an average observed formation time of 0.6 Myr, cor- 
responding to a mean disk lifetime of about 3 Myr, at 
the high end of the observationally inferred values. Im- 
provements in the accuracy of the measured protostellar 
lifetimes, acceleration times and luminosities will enable 
more stringent tests of the models. 

4.5. The Protostellar Luminosity Function (PLF) 

In this section we present the PLF for each model and 
define four dimensionless parameters to characterize the 
PLF shape. In Figures 5 and 6, we plot the predicted 
PLFs and overlay the PLF of the observed distribution 
of protostars. We exclude sources with L < 0.05 Lq 
(comparable to the limit of the observations after cor- 
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TABLE 2 

Median Luminosity (Lmed(i©) with (t/)^^^ = 0.44 Myr; 

-^med,obs — — U.4 -^0 ■ 



Non- Accelerating Accelerating" 
Model Untapered Tapered'' Untapered Tapered*' 



Isothermal Sphere (IS) 3.12 

2C Turbulent Core (2CTC)'= 1.76 

Turbulent Core (TC) 1.44 

2C Competitive Accretion (2CCA) 0.95 

Competitive Accretion (CA) 0.88 

= 1 Myr 
''n = 1 
'^Tl^ = 3.6. 

rectioii for extinction) from the model distributions and 
then renormalize the PLF to unity. 

In Table 3 we give the standard deviation of log lumi- 
nosity, the ratio of the median luminosity to the mean, 
the ratio of the maximum luminosity to the mean, and 
the fraction of very low luminosity objects (VELLOs; see 
Section 4.5.4) for each model, re-normalized such that 
the means are equal to the observed mean luminosity. 

4.5.1. Standard Deviation of Log Luminosity, (j(logL) 

We find that the IS model has the smallest cr(logL) 
in all cases and is too narrow with respect to the data. 
In contrast, the TC and CA models have larger values 
of (7(logL) and encompass the extent of the data. The 
2CTC and 2CCA models present a promising compro- 
mise and their cr(logi) are slightly too low. Allowing 
the accretion rate of the protostars to taper off decreases 
the standard deviation of the distributions in all but the 
non-accclcrating IS case. This exception occurs because 
tapering introduces a spread in the accretion rate. 

Allowing for time variations of a factor of two, i.e., non- 
FU-Ori fluctuations in the accretion rates, would also 
increase the standard deviation but have little effect on 
{L) or L-n-icd/ {L) ■ To take such low-level variability into 
account, we add 0.3 dcx in quadrature to the measured 
standard deviations. As shown in Table 3, this permits 
agreement between observations and both the 2CTC and 
2CCA models. The TC and CA models remain within 
error, while the IS models continue to be inconsistent. 

In Figure 8, we plot (T(logL) as a function of the most 
massive forming star. We plot the tabulated values for 
the case of m„ = 3 Mq and the observed standard devia- 
tion {(j{\ogL) = 0.7^Q^ dex) in inset plots. As illustrated 
by both Figures 8 and 5, the standard deviations of the 
four models are very distinct. The TC and CA models 
grow closer together for larger rn„, but remain fairly well 
separated at m„ = 3 Mq. The standard deviations are 
also relatively insensitive to the upper stellar mass in the 
cluster. Tapering reduces the width of the distribution 
for all m„ for these models. 

4.5.2. Median to Mean Ratio 

The ratio of the median to the mean, L^q^/ {L), vir- 
tually eliminates the dependence of the results on {tf), 
/epi, and /acc- Table 3 gives the values of the ratios for 

each of the accretion models, where the observed ratio 
is 0.3^t^. In all cases, the IS model can be ruled out.^ 

^ Dunham et al. (2010) find that the IS model as characterized by 
constant accretion onto a disk may still be consistent provided that 



3.19 


3.86 


3.22 


1.85 


1.41 


2.44 


1.52 


1.91 


2.22 


1.75 


1.52 


2.76 


1.25 


1.31 


1.69 



The non-accelerating 2CTC and the tapered 2CCA mod- 
els are outside the observational error bars. 

4.5.3. Maximum Luminosity, Lmax 

To characterize the maximum luminosity of the distri- 
butions, we define: 

ima. = y " *(iacc) dhlL = 1/N, , (43) 

where >I'(Lacc) is the PLF derived using the accretion lu- 
minosity. For parity between the observations and mod- 
els, we compare with the ratio of I/max to (Lobs)- The 
most luminous observed source has a bolometric lumi- 
nosity of 76 Lq, which is likely an upper limit and may 
be over estimated by a factor of ^ two. Consequently, we 
adopt 54^^g Lq in comparing with the models. Impor- 
tantly, this source is a Class object, which suggests that 
the luminosity chiefly arises from accretion. Thus, we de- 
rive Lniax from ^'(Lacc), which assumes that the interior 
luminosity is small compared to the accretion luminosity. 
We note that the values of Lmax in Table 3 are based on 
the actual model dispersion, a, not the enhanced value 
(Teff that allows for factor of 2 variability; allowance for 
such variability would increase Z/max somewhat. 

According to Table 3, the untapered CA models have 
the highest imax/(-Zjobs), followed by the . 2CCA and 
2CTC cases and the TC ones. The observed maximum 
luminosity is ~ 2 — 3 times higher than for the tapered 
PLFs, but it is within error of the untapered 2CTC, 
2CCA, TC and CA cases. Without a correction for the 
model mean luminosities, the difference between the ob- 
served and model Lmax would be a factor of 5-6. It is un- 
clear whether this discrepancy is an artifact of the larger 
error of the higher luminosity measurements, variable ac- 
cretion, or fluctuations in the star formation rate. 

The untapered cases, even though directly consistent 
with Lmax/ {L), are discrepant with the observations in a 
different way. For untapered accretion, in which the ac- 
cretion rates increase monotonically, the maximum lumi- 
nosities are achieved only when m reaches its final value. 
rrif. Consequently, in order for the untapered cases to 
be consistent with observations either the Class I phase 
must have higher observed luminosities than the Class 
phase or the Class phase must last much longer is 
currently inferred from observations. The tapered IS, 
2CTC, and TC models, though lower as shown by Ta- 
ble 3, have peak luminosities that occur midway through 

the accretion from the disk onto the star is completely episodic. 
Here, we rule out the IS model where accretion from the disk onto 
the star is smooth. 
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TABLE 3 
Model Results for (L) = (Lobs) 



Model 




(T(log L) 


(Teff (log LY 




Lmax (Class 0)/(L) 


/VBLLO'' 


(t/>-(t/>ob.(Myr)= 


Observed 








3+°-2 


lOl^ 


0.2 ±0.1 


0.44 ±0.22 


Isothermal Sphere (IS) 


Non-Accclcrating 
Untapcrcd 
Tapered'* 

Accelerating'^ 
Untapcrcd 
Tapered 


0.44 
0.45 

0.43 
0.39 


0.55 
0.56 

0.54 
0.51 


0.96 
0.83 

0.86 
0.71 


5.81 
2.97 

2.90 
3.67 


0.08 
0.08 

0.07 

0.07 


0.62 
0.72 

0.85, 0.37 
0.85, 0.37 


Two-Component 
Turbulent Core (2CTC)f 


Non-Accclcrating 

Untapcrcd 

Tapered 
Accelerating 

Untapered 

Tapered 


0.55 
0.41 

0.56 
0.42 


0.64 
0.53 

0.65 
0.53 


0.60 
0.64 

0.26 
0.45 


11.80 
5.23 

10.57 
6.74 


0.09 
0.08 

0.15 
0.07 


0.56 
0.60 

0.99, 0.43 
1.03, 0.45 


Turbulent Core (TC) 


W on- Accelerating 

Untapered 

Tapered 
Accelerating 

Untapered 

Tapered 


0.77 
0.69 

0.76 
0.71 


0.84 
0.76 

0.83 
0.78 


0.42 
0.51 

0.32 
0.53 


9.46 
4.91 

9.78 
6.06 


0.18 
0.21 

0.21 
0.12 


0.69 
0.58 

1.08, 0.47 
0.84, 0.37 


Two-Component Competitive 
Accretion (2CCA)8 


Ixlon— A cr'fA (^va + in cr 

Untapered 
Tapered 
Accelerating 
Untapered 
Tapered 


0.53 
0.49 

0.55 
0.48 


0.61 
0.57 

0.63 
0.57 


0.30 
0.66 

0.30 
0.67 


9.46 
4.04 

12.19 
4.02 


0.09 
0.06 

0.11 
0.05 


0.55 
0.51 

0.94,0.41 
0.78,0.34 


Competitive Accretion (CA) 


Non-Accclcrating 

Untapcrcd 

Tapered 
Accelerating 

Untapered 

Tapered 


0.81 
0.75 

0.80 
0.80 


0.94 
0.89 

0.93 
0.93 


0.21 
0.44 

0.18 
0.43 


12.78 
5.99 

13.41 
6.79 


0.25 
0.22 

0.27 
0.17 


0.81 
0.62 

1.33, 0.58 
0.89, 0.39 



^Effective standard deviation, which assumes a fax;tor of two time-variability in the observed luminosities; 0.3 dex is added in quadrature to the 
model values given for iT(log L). 
''See equation 44. 

'^For accelerating models, both the actual mean star-formation time, (t/), and the one that would be inferred by comparing the number of 
protostars with the number of Class II sources, {tf)^^^, are given. These times are equal for non-ax;celerating models. 
•^n = 1 
'=T = 1 Myr 
'TlrH = 3.6 
^T^m.CA = 3.2 



the formation time. In contrast, the untapered 2CCA 
and CA models, which agree better with the observa- 
tions, achieve those luminosities only towards the end of 
accretion, and thus do not appear to be consistent with 
observation. 

4.5.4. Very Low-Luminosity Object Fraction, /vello 

There is observational evidence in support of a signifi- 
cant population of low luminosity protostars. For exam- 
ple, Dunham et al. (2008) found that ^ 30% of protostars 
with L < 1.0 Lq have luminosities below 0.1 Lq, a sam- 
ple commonly referred to as very low-luminosity objects 
(VELLOs). In our sample, which covers the same re- 
gions as Dunham et al. (2008) but has been corrected for 
dust extinction and more carefully pruned to eliminate 
non-protostars (e.g., background galaxies), only ^ 20% 
of protostars can be considered VELLOs. Note that ap- 
plying an exinction correction to the data increases the 
median luminosity by 40% so we define the VELLO frac- 
tion as 

The observational sample contains one source with an 
extinction-corrected luminosity of 0.035 Lq, but it is 



likely incomplete at luminosities below 0.05 Lq; we there- 
fore set Lmin = O.O5L0. Without correcting for ex- 
tinction, ~ 20% of the observed sources would have 
L < O.ILq. 

In Table 3 we give /vello for each of the models. The 
IS models, the tapered 2CTC, and the tapered 2CCA 
models are inconsistent with the data, whereas the TC 
and CA models are consistent with the observed values 
in all cases. The CA, and to a lesser extent, the TC mod- 
els also predict a substantial mimber of VELLOs below 
the luminosity completeness limit. Future observations 
extending below this limit are necessary to confirm or 
rule out these models. It should be noted that in all 
prescriptions, VELLOs are produced not by quiescent 
periods prececding episodic accretion events, but by the 
small accretion rates associated with protostars of very 
low mass. Current observations cannot discriminate be- 
tween low-luminosity sources in a quiescent phase and 
those that are simply low-mass objects with low accre- 
tion, although it is likely that at least some of the ob- 
served protostars fall in the former category. It is there- 
fore possible that the low VELLO fractions of some of 
the models may be consistent with the actual fraction of 
the subset of very low-luminosity protostars that are not 
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Fig. 3. — The mean total protostoUar luminosity versus the upper 
protostellar mass, m„, for each of the models with (tj)= 0.44 Myr. 
For the tapered cases, n = 1, and for the accelerating cases, r = 1 
Myr. The inset shows the mean observed luminosity from Evans 
et al. (2009) with error bars representing the uncertainty in the 
measurement and rrtu. The diamonds display the values in Table 
1. These values assume Lmin = 0.05 Lq and (tf)^^^ = 0.44 Myr, 
which corresponds to {tf) = 0.24 Myr for the accelerating cases. 
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Fig. 4. — Protostellar lifetime estimated using the observed mean 
luminosity from the Evans et al. (2009) data as a function of the 
mean luminosity from the models. The error bars on the model 
lifetimes derive from the observal luminosity uncertainty. The two 
observational results with uncertainty are shown by the thick set 
of solid error bars. 
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Fig. 5. — The PLF for rrtu = SM© for untapered, non- 
accelerating star formation (top) and untapered, accelerating star 
formation with t = 1 Myr (bottom). The observed PLF (Evans 
et al. 2009) is plotted for comparison. Note that the PLF shape 
is derived assuming that the accretion luminosity dominates the 
total. 

undergoing episodic accretion. 

4.5.5. Constant Radius PLF 

In the Appendix we derive the PLF for each model 
(except 2CCA) assuming only accretion luminosity and 
adopting a constant protostellar radius, r. Figure 7 
shows these constant-radius curves in the untapered, 
non-accelerating case together with the PLFs for which 



The Protostellar Luminosity Function 



13 



0.6 



0.4 



0.2 - 



0.0 



All protestors 
IS 

2CTC 

TC 

CA 



Topered, 

Non-occeleroting i — y 




0.1 



1.0 10.0 

L 



100.0 



0.6 



0.4 



0.2 



0.0 



All protostors 
IS 

2CTC 

TC 

CA 



Topered, 
Acceleroting 



n 




0.1 



1.0 10.0 
L 



100.0 



Fig. 6. — The PLF for rrtu = 3A/0 with tapered accretion rates. 
The observed PLF (Evans et al. 2009) is plotted for comparison. 
Note that the PLF shape is derived assuming that the accretion 
luminosity dominates the total. 
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Fig. 7. — The PLF for mu = 3Mq for the fiducial cases, where 
r is a constant and where r is a polynomial fit to f(m). In each 
model, the curves are normalized to the same mean luminosity, (L) 
(see §4.2). 

r is a polynomial fit to the sub-grid protostellar evolu- 
tion model (as in Figure 5). Figure 7 illustrates that 
the curve width and maximum luminosity are reduced 
when the radius is allowed to depend upon mass. The 
constant radius curves have standard deviations that are 
^ 50% larger than when the radius is allowed to vary, a 
change that, for example, makes the shapes of the con- 
stant radius IS curve and varying radius 2CTC curves 
similar. This indicates that the stellar evolution model 
does impact the PLF shape and the details of the com- 
parison. In particular, stellar evolutionary models with 
a narrower range of r will have broader luminosity distri- 
butions. Since our stellar evolution model depends upon 



the instantaneous accretion rate and not just the accre- 
tion timescale, a star of the same final mass will have a 
different amount of deuterium remaining at t = tf for 
the different accretion histories. The models fall in the 
order IS, 2CTC, 2CCA, TC, and CA from most evolved 
to least evolved, based upon the amount of deuterium 
burned at a given final mass. However, we note that 
the model PLF shapes are mainly distinguished by their 
characteristic accretion models rather than differences in 
their stellar evolution. 

5. DISCUSSION 

A number of uncertainties enter into our comparisons. 
Foremost, the observational data are difficult to obtain 
and corrections due to obscuration along the line of sight 
contribute significant error. There is also uncertainty in 
several of the parameters we have adopted for the models 
and in the way these parameters are implemented. 

The parameter /acc remains one of the most uncertain 
in our estimation. It actually contains two very differ- 
ent effects: non-radiative loss of accretion energy from 
the disk (e.g., Ostriker & Shu 1995) and advection of ac- 
cretion energy into the stellar interior (Hartmann et al. 
1997). The latter effect reduces the accretion luminos- 
ity by a factor (1 — a), where a is the fraction of the 
accretion energy advected into the stellar interior. Most 
authors agree that this effect is small (Hartmann et al. 
1997; Baraffe et al. 2009), and it was not included in the 
work of Stabler (1988). In a thorough study of the struc- 
ture of protostellar accretion shocks, Commergon et al. 
(2011) have shown that a is indeed extremely small pro- 
vided the accretion fiow is optically thin to optical radia- 
tion, so we neglect it here. The accretion shock does heat 
the surface layers of the protostar, however, and if one 
assumes that this heating is negligible because the ac- 
cretion is localized onto a small fraction of the protostel- 
lar surface, the resulting protostars can be significantly 
more compact than when the surface of the protostar is 
heated by the accretion shock (Hartmann et al. 1997). 
Such models appear to be inconsistent with observation 
(Baraffe et al. 2009; Hosokawa et al. 2011). 

The fraction of mass accreted during bursts, /epi, is 
also uncertain. Most bursts are observed to occur in the 
Class I stage, although Class II objects also experience 
FU-Ori type events (e.g.. Miller et al. 2010). There is 
some suggestion that bursts increase with age and that 
Class objects experience smoother accretion (Vorobyov 
& Basu 2005; Zhu et al. 2009). There is also evidence 
that isolated stars experience more frequent outbursts 
than those in clusters (Greene et al. 2008), a puzzle that 
highlights a deficit in our understanding of protostellar 
accretion processes. 

If protostars undergo periodic outbursts, then they 
must exist in a quiescent, low-luminosity stage between 
bursts. Indeed, some of these quiescent sources may ap- 
pear as VELLOs. The median luminosity of the FU 
Ori sources cataloged by Sandell & Weintraub (2001) is 
250 Lq, so protostars that increase in luminosity by more 
than 8 magnitudes would have been in a VELLO state 
prior to outburst. It is not known how many of the known 
FU Ori sources had pre-burst luminosities below 0.14 Lq, 
but several of the best known ones did have pre-burst lu- 
minosities above this limit (Hartmann & Kenyon 1996) 
and would not have been classified as VELLOs. If the 
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Fig. 8. — The standard deviation of the log of the protostellar 
luminosity as a function of the upper protostellar mass, m„, for 
the models with (tf) = 0.44 Myr. For the tapered and accelerating 
cases n = 1 and t = 1 Myr, respectively. In the inset, the standard 
deviation of the observed luminosities (Evans et al. 2009) is plotted 
with error bars to represent the uncertainty in the measurement 
and mu- The diamonds display the values from Table 3, which 
assume Lniin = 0.05 Lq and {tf)^y^^ = 0.44. 



bursting sources were in a low-luminosity state between 
bursts, they would populate the lower luminosity region 
of the observed PLF and thus increase the standard devi- 
ation, lower the ratio of the median to the mean, and, for 
sufficiently faint sources, increase /vello of the observed 
population relative to our models, which do not explicitly 
include outbursts. If the inter-burst accretion luminosity 
is negligible and the nuclear luminosity is approximated 
by the ZAMS value, then non-accreting protostars with 
masses up to about 0.7 Mq could potentially be classi- 
fied as VELLOs. The observed value of /vello is thus 
an upper limit on the number of sources in that lumi- 
nosity range in our models. As a result, models like the 
tapered 2CTC and 2CCA models, which fall below the 
observed /vello, are promising candidate models, while 
the untapered TC and CA accretion models may well 
over-predict /vello and are not as promising. 

The ratio of the median and mean luminosities can be 
significantly reduced by the degree of burstiness in the 
accretion rate of the individual protostars. The same ef- 
fect can be achieved in models such as the CA and TC 
models, which achieve a small ratio by virtue of a rela- 
tively long phase of low-luminosity accretion. However, 
as noted above, these models are somewhat artificial, 
the former because protostars initially accrete some mass 
from a local reservoir (which leads to the 2CCA model) 
and the latter because turbulent motions do not exceed 
thermal motions in such low-mass cores (which leads to 
the 2CTC model). Nonetheless, this illustrates that it 
is possible to fit the observed ratio without episodic ac- 
cretion, and thus, that it is observationally difficult to 
disentangle the infiuence of variability from an underly- 
ing mean accretion trend shaping the luminosity distri- 
bution. 

Although only a small number of sources have been 
observed to undergo large, FU Ori outbursts, many pro- 
tostellar objects have been observed to undergo variabil- 
ity in luminosity over a factor of two on timescales of 
months to years (Pech et al. 2010; Covey et al. 2010). 
Our correction to the dispersion in log L discussed in 
§4.5.1 reflects this, but is very approximate. In fact, we 
find that some low-level time-variability is needed in or- 
der for the 2CTC and 2CCA cases to be consistent with 
observations. 

Another source of dispersion in the PLF could stem 
from variation in the initial conditions for star formation. 
For example, temperatures in the low-mass star-forming 
regions we study here are in the range 10 — 20 K (McKee 
& Offner 2010a), which leads to a ~ factor 3 range of 
accretion luminosity in the IS model. This corresponds 
to a dispersion of ~ 0.15 dex, which has a negligible effect 
when added in quadrature to the dispersion of 0.3 dex for 
the assumed temporal variability. When samples with 
a larger range of initial conditions are considered, it is 
possible that this additional source of dispersion would 
have to be included. 

Our models do not directly take into account other 
modes of star formation such as fragmentation within 
accretion disks (Bate 2009a; Stamatellos & Whitworth 
2009; Kratter et al. 2010). Division of gas accreting onto 
a shared disk between close companions is a complicated 
process, which could potentially alter the the accretion 
dependence on m and mf that we assume. However, 
the formulation of the models does not necessarily ex- 
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elude disk fragmentation since we specify the accretion 
prescription rather than the protostellar origin. We also 
do not expect disk fragmentation to be common in this 
observational sample, since heating due to radiative feed- 
back significantly stabilizes low-mass disks (Offner et al. 
2009; Bate 2009b). In the TC and CA senarios, core 
and filament fragmentation may produce wide binary 
companions (Offner et al. 2010), but accretion rates for 
this mode of origin should follow the expected accretion 
trends. We have not taken binarity into account in our 
analysis, but its effects are small compared to the large 
differences among the different accretion models; for ex- 
ample, the difference between Chabrier's (2005) IMF for 
individual stars and that for stellar systems is only a 
factor of 1.25 in the peak mass. 

In addition to the accretion history and to uncertain- 
ties in parameters, inclination effects may also contribute 
to the standard deviation of the observed PLF. For ex- 
ample, the extinction correction and therefore the bolo- 
metric luminosity of a protostar with a disk that is ob- 
served edge-on will be underestimated. Dunham et al. 
(2010) found that adding inclination effects to their mod- 
els broadened the bolometric luminosity distribution by 
< 20%, so that orientation effects alone were not able to 
fit the data. Moreover, including inclination effects had 
only a small effect on the high luminosities, at most in- 
creasing the maximum by 25%. Consequently, we expect 
geometry to have a minor effect on the shape of the PLF. 

Adopting an accelerating star formation rate is one 
possible resolution of the apparent discrepancy between 
the observed star formation timescale of 0.44 Myr, which 
is based on the assumption of a constant star-formation 
rate (Evans et al. 2009), and the mean luminosity, which 
suggests a star-formation timescale of about 0.3 Myr 
(Sec. 4.4). We have adopted r = 1 Myr as the accel- 
eration time, which is the inferred acceleration time for 
Ophiuchus and is likely to be a lower bound on the accel- 
eration time for the other regions (Palla & Stabler 2000; 
McKee & Offner 2010b). (More recent Ophiuchus data 
suggest that the region has a deficit of Class objects 
and therefore a decelerating rate of star formation (Evans 
et al. 2009).) Insofar as 1 Myr is a lower bound on the ac- 
celeration time, it gives the maximum difference between 
(tf) and (i/)^!^^, and thus the maximum effect of acceler- 
ation on protostellar luminosities. IC348, a sub-region of 
Perseus, has an inferred accleration time of 2 Myr (Palla 
& Stabler 2000; McKee k Offner 2010b), but the accel- 
erations for the entire Perseus region and for Serpens are 
unknown. Given the small statistical sample size, the 
appreciable uncertainties in the timescale estimates, and 
our simplified acceleration model, these results should be 
interpreted as being consistent with, but not demonstrat- 
ing that, acceleration resolves the apparent discrepancy 
between the observed star-formation time and the mean 
luminosity by accelerating star formation. 

How do the different accretion models compare with 
the data? Independent of parameter uncertainties, the IS 
models appear to be inconsistent with observations. Of 
the five metrics we have adopted — (Toff(logL), £mcd/ {L), 
/vELLO, imax/(i), and (tf)^^^— Only the last is within 
the estimated errors. The ratio of the median to mean 
luminosity is more than 2 standard deviations from the 
observed value. (However, the error estimates here are 
quite uncertain, and future data should substantially im- 



prove them.) Because protostars accreting according to 
the TC and CA prescriptions spend significant time at 
low masses, where IS-like accretion is likely to occur, we 
believe that the 2CTC and 2CCA models are more real- 
istic physically, even though in some cases the TC and 
CA models agree quite well with the data. For example, 
in the CA case low- mass stars spend 1/3 of their accre- 
tion time achieving a mass of 0.1 Mq. It is unlikely that 
protostars spend this much time at accretion rates sig- 
nificantly less than the IS value. The 2CTC and 2CCA 
models are thus likely to be a better representation of the 
TC and CA models, since the constant accretion compo- 
nent dominates the initial protostellar accretion rate. 

Likewise, with the exception of the CA models, which 
are assumed to have accretion terminated over a short 
time by protostellar feedback, tapered models are likely 
to be more realistic since accretion from a core; dciclincis 
gradually after the expansion wave reaches the surface 
of the core (e.g., McLaughhn & Pudritz 1997). Further- 
more, untapered models achieve their maximum lumi- 
nosity at the end of the protostellar stage, whereas ob- 
servations show that Class I sources are not noticeably 
more luminous than Class ones; in addition, the accre- 
tion rate for the Class sources is believed to be larger 
than for the Class I sources, which is consistent with ta- 
pering. In our comparison we find that in all cases the 
tapered models tend to underestimate the maximum lu- 
minosity and the standard deviation of the distribution. 
This may suggest that our parametrization of tapering 
or our specific choice of n = 1 is not correct rather than 
ruling out tapered protostellar accretion rates. It must 
also be borne in mind that we have not included the ef- 
fect of variability on Lmaxj so the tabulated values are 
lower bounds on the true values. 

As a check on our tapering model, we derive tapered 
PLFs for the IS and TC cases assuming an exponential 
functional form for the accretion rate: 

m = mo{m,mj)cxp{-2t/tf) {t<tf). (45) 

For the IS case, this is similar to the model of Myers 
et al. (1998), except that in their model m approaches 
m,f as t ^ OG. We modify their model so that m. = m,f 
at t = tj\ at which point the accretion rate has declined 
exponentially by a factor of e~^. For the IS case, we 
find that cr(Log L), (L), L^^d/ (L) and L^^ax/ (L) change 
by only a few percent relative to the linearly tapered 
case. The parameter /voiio decreases by 50%, mainly be- 
cause the exponentially tapered accretion rate does not 
go to ~ as the linearly tapered case; does. Thus, adopt- 
ing a different tapering model does not alter the fun- 
damental disagreement of the IS case with the observa- 
tions. For the exponentially tapered TC case, a(Log L) 
and £max/(-^) change by a few percent, while (tf) and 
Lmcd/{L) increase by ~ 20 and ^ 25%, respectively, and 
/veiio decreases by 40%. Although these changes are more 
significant, the dimensionless metrics continue to agree 
with the observations as before, except that L-a^^i/ {L) be- 
comes slightly too high. While it may be possible to for- 
mulate accretion tapering that improves the agreement 
with the observations, it seems unlikely that a different 
functional form would modify our conclusion that con- 
stant accretion time models agree better with the data 
than constant accretion rate models. 
The comparison with observations could be strength- 
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Fig. 9. — The PLF for the untapered, non-accelerating star for- 
mation 2CCA model with ca = 2.0, 3.6 and the corresponding 
IS, CA, and observational PLFs for comparison. 

ened by deriving PLFs for the Class and Class I pop- 
ulations separately. In theory, the division between the 
two classes occurs when the protostellar mass exceeds 
the envelope mass (Andre et al. 2000; Crapsi et al. 2008). 
For the IS and TC models, it is straightforward to de- 
fine Class and Class I PLFs as the distributions of lu- 
minosities for which m < \mf and i^mj < m < nif, 
respectively. Although such distributions can also be 
constructed for the CA case using our approximate CA 
model, one characteristic of competitive accretion is the 
lack of correspondance between envelope mass and final 
stellar mass (Smith et al. 2009). Even using this simple 
definition, comparing the Class PLFs quantitatively with 
the observations is challenging. Since the protostellar 
mass is not measurable during the embedded phase and 
since inclination effects can significantly confuse the clas- 
sification, the two populations are difficult to distinguish 
observationally. Because the Class I lifetime is about 
three times the Class lifetime (Evans et al. 2009), mod- 
els that do not have a higher rate of accretion as Class 
sources could be excluded. If the fraction of a protostel- 
lar lifetime spent as a Class source were to be less than 
that found by Evans et al. (2009), the constraint on ac- 
cretion and early luminosities would become more strin- 
gent. Qualitatively, we expect all the untapered models 
to fail such a test, since accretion, and hence luminosity, 
in such models rises monotonically with age. 

6. CONCLUSIONS 

We have analyzed the Protostellar Luminosity Func- 
tion (PLF), which is the present-day luminosity distri- 
bution of a cluster of protostars, for several different the- 
ories of star formation and compared the results with 
observation. In our derivation we have assumed that 
the protostellar masses evolve smoothly onto a truncated 
Chabrier (2005) IMF, that the accretion rates are a con- 
tinuous function of the instantaneous and final proto- 
stellar masses, and that the star formation rate is either 
constant or accelerating in time. We also assume that 
over most of the formation time, the accretion rate onto 
the protostar tracks the accretion rate from the ambient 
molecular core onto the protostar-disk system. Episodic 
accretion, such as that occurring in FU Ori outbursts, 
violates this assumption, but the available observational 



data indicate that this is not a dominant effect. 

The PLF depends explicitly upon the mean formation 
time of the stars, (t/), and the maximum stellar mass 
produced, m„. For the low- mass star- forming regions 
we have compared with, m„ appears to be about 3Mq. 
We consider three main accretion prescriptions corre- 
sponding to standard models of star formation: Isother- 
mal Sphere (IS, constant accretion rate). Turbulent Core 
(TC), and Competitive Accretion (CA, constant accre- 
tion time). We note that prior to the development of 
either the CA or TC models, Kenyon et al. (1990) con- 
sidered both constant accretion rate and constant accre- 
tion time in the paper that introduced the luminosity 
problem. We also consider two hybrid models: the Two- 
Component Turbulent Core model (2CTC, a compromise 
between the IS and TC models) and the Two-Component 
Competitive Accretion model (2CCA, a combination of 
the IS and CA accretion prescriptions; this model was 
not considered in Paper I). We explore two variations on 
these models: a case in which the accretion rate smoothly 
tapers to zero and a case in which the star formation rate 
accelerates with a characteristic timescale, r. 

The CA model used here is an approximate analytic 
representation of the competitive accretion model devel- 
oped by Bonnell et al. (1997), which begins with proto- 
stellar seeds that are produced by a process similar to 
that in the IS case. As a result, we believe the 2CCA 
model is a better approximation to their work than the 
CA model. The TC model was specifically formulated for 
high-mass stars, which we do not focus on here. For the 
case of low-mass stars, McKee & Tan (2003) proposed the 
inclusion of an IS stage, which suggests that the 2CTC 
model is the best representation of their model in this 
context. We note that this model has some similarities 
to the TNT model of Fuller & Myers (1993). 

We compare our models to protostellar luminosities ob- 
served in local low-mass star forming regions (Evans et 
al. 2009, Enoch et al. 2009). The classical luminosity 
problem is that observed protostars appear to have lu- 
minosities significantly lower than expected theoretically 
(Kenyon et al. 1990). The extinction-corrected sample 
that we adopt from Evans et al. (2009) has a mean lumi- 
nosity of 5.3l^'g Lq, which is more than a factor of two 
larger than the earlier Enoch et al. (2009) results that 
did not account for extinction. This alone signficantly 
ameliorates the luminosity problem. In comparing the 
models and the data, we first used the mean luminos- 
ity and the median luminosity. In all permutations of 
the parameters, the models require that the average pro- 
tostellar lifetime be 0.3 ± 0.1 Myr, somewhat less than 
> 0.4 Myr measured by Evans et al. (2009). That is, 
the model luminosities are too low if the mean lifetime 
inferred by Evans et al. (2009) is assumed. Thus, with 
a star formation time of 0.3 Myr, allowance for non- 
radiative energy loss in winds (/acc — 0.25) and a modest 
amount of episodic accretion (/cpi — 0.25) is sufficient to 
lower the mean protostellar luminosity so that there is no 
longer a "luminosity problem," in low-mass star forma- 
tion. We note that this resolution of the luminosity prob- 
lem is quite consistent with the suggestions of Kenyon 
et al. (1990), who first pointed out the existence of the 
problem: among the solutions they proposed were that 
the formation time was longer than the (1 — 2) x 10^ yr 
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indicated by their data and that some of the accretion 
was episodic. With a star-formation time of 0.3 Myr, the 
mean accretion rate for a star of mass 0.5Mq (the mean 
mass oftheChabrier 2005 IMF) ism ~ 2xlO-^M0 yr'^ 

The discrepancy between our estimate of the mean star 
formation time and that of Evans et al. (2009), while not 
large, could be due to a number of factors: the disk life- 
time could be shorter than they assumed, the number of 
Class I sources could less than they assumed, fepi or /acc 
could be larger than we assumed, and/or the star forma- 
tion rate could be accelerating as suggested by Follow- 
up studies of dense gas tracers, such as HCO^, have 
found that as many as half of previously identified Class 
I objects arc not actually embedded van Kempen et al. 
(2009); Heiderman et al. (2010). Although Evans et al. 
(2009) exclude sources embedded in less than 0.1 Mq of 
gas mass when deriving the lifetime, some sources may 
nonetheless be misclassified older, heavily obscured ob- 
jects. Also, while our value of /epi is calculated using all 
the known bursting sources, deeply embedded FU-Ori 
type sources may be missing from the sample. By in- 
creasing the ratio of protostars to Class II sources, accel- 
erating star formation produces a longer observationally 
inferred star formation time than the actual star forma- 
tion time of the sample. Consequently, the accelerating 
models have better consistency between the observed star 
formation time and mean and median luminosities. Al- 
though we do expect the star formation rate to vary with 
time, Palla & Stabler (2000) found evidence for an accel- 
eration time as short as our adopted value, r = 1 Myr, 
for only one the observed regions in the protostellar sam- 
ple (Ophiuchus), and this evidence is contradicted by an 
apparent deficit of Class O sources in the region. 

We then compared four dimensionless quantities that 
characterize the shape of the PLF: (1) the standard de- 
viation of log L (including an allowance for source vari- 
ability at the factor of two level); (2) the ratio of the 
median to the mean luminosity; (3) the ratio of the max- 
imum to the mean luminosity; and (4) the fraction of 
very low-luminosity objects, /vellOi defined as the ra- 
tio of the number of sources with extinction-corrected 
luminosities between O.OSL© and 0.14_Lq to the num- 
ber between 0.05Lq and 1.4^©. The first three of these 
are the most strongly discriminating since they arc inde- 
pendent of /epi, /acc and the mean protostellar lifetime, 
(tf); the fourth quantity, /vellO) is only weakly depen- 
dent on these factors. We also compared the value of 
the star-formation time required by the models to get 
the observed mean protostellar luminosity with the star- 
formation time of = 0.44 ± 0.22 Myr inferred by 
Evans et al. (2009) from the ratio of the number of pro- 
tostars to the number of Class II sources, which were 
assumed to have a 2 Myr lifetime. Although protostars 
with different accretion histories have slightly different 
stellar evolutionary states at the end of accretion, we find 
that differences between the PLF shapes are driven by 
the different accretion histories, not the different evolu- 
tionary states. We find that the IS model is a poor fit to 
the data in all cases, mainly due to the strongly peaked 
nature of its PLF profile. The model results for the four 
dimensionless parameters are outside the error bars of 
the data regardless of whether the model is tapered or 
untapered, accelerating or non-accelerating. The one pa- 



rameter the IS model agrees with is the observed star- 
formation time, after renormalizing the accretion rate so 
that the mean luminosity agrees with observation. 

The CA model is in best agreement with the data; only 
the observed star-formation time, {tf)^^^, lies outside 
the error bars, and this is only for the non-accelerating, 
untapered case. The CA model predicts a relatively 
large number of VELLOs below the observational limit 
of 0.05^0, which will provide a strong test of the model 
in the future. However, as discussed above, the 2CCA 
model, which has an initial phase in which the star ac- 
cretes more rapidly, is closer to the actual competitive 
accretion model. Furthermore, one of the assumptions of 
the competitive accretion model is that the gas is cleared 
out relatively quickly toward the end of the accretion 
phase, so the untapered version of the 2CCA model is 
closest to the actual competitive accretion model. How- 
ever, this model (as well as the untapered CA model) 
predicts that the maximum himinosity is achieved at late 
times, not during the Class stage. In general, the length 
of the Class lifetime and the similarity of Class and 
Class I luminosities suggests that models that do not 
accrete a significant fraction of mass during the earliest 
times are inconsistent. This includes all the untapered 
models, which acheive their maximum accretion rate, 
and hence maximum luminosity, at the end of the proto- 
stellar lifetime, in what would be the late Class I stage. 
Otherwise, both the accelerating and non-accelerating 
untapered versions of the 2CCA model agree well with 
the data, although the dispersion is very slightly below 
the observed value for the non-accelerating case. The ta- 
pered version of the 2CCA model does not compare well 
with the data. 

The TC models are in good agreement with the data, 
with the exception of the non-accelerating, tapered case, 
which has a maximum-to-mean luminosity ratio that is 
marginally too low. The observed star-formation time, 
(*/)obs' marginally too high for the non-accelerating, 
untapered case. As for the CA model, however, the 
2CTC model, which is similar to the IS model at low 
masses, is more realistic. The untapered version of this 
model agrees well with the data, except that the median 
luminosity is somewhat high for the non-accelerating 
case. However, the tapered version of the 2CTC model 
is the best representation of the model, since the model 
is essentially a turbulent version of the IS model. The 
non-accelerating version of this model marginally agrees 
with all the data except for Lmax/(i); the accelerating 
version underprcdicts the number of VELLOs (although 
as remarked above, this may not be a problem if some of 
the observed VELLOs are episodic sources in a quiescent 
phase). 

We conclude that models that tend towards a constant 
accretion time and thus produce a greater spread in lu- 
minosities (like the CA and TC models), rather than 
models that have a constant accretion rate (such as the 
IS model) are in better agreement with the data on the 
PLF. Ultimately, agreement between a model and the ob- 
served PLF is necessary but not sufficient. Models must 
also reproduce a number of other observed features of 
protostars and the regions from which they form, includ- 
ing core properties (e.g. Offner et al. 2008; Kirk et al. 
2009), the approximate agreement between the luminosi- 
ties of Class and Class I sources, and the existence of 
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large disks around Class I sources. The CA model, whieh 
exhibits good agreement with the PLF, does not do well 
with any of these additional features; in particular, it 
does not produce well-defined cores for individual proto- 
stars (Enoch et al. 2008). The IS and TC models natu- 
rally have cores, and the tapered models can give com- 
parable luminosities for the Class and Class I stages. 
Non-magnetic IS and TC models are predicted to have 
large disks, but the existence of such disks in the pres- 
ence of magnetic fields is a topic of active investigation 
(e.g., Mellon & Li 2009; Ciardi & Hennebelle 2010). 

In this paper, we have developed the PLF as a tool 
for confronting star formation theories with observa- 
tion. Ongoing observational efforts should permit a sig- 
nificant improvement in the comparison between the- 
ory and observation by providing larger samples, which 
would reduce statistical fluctuations, and more accu- 
rate extinction-corrected luminosities, which would re- 
duce the current factor of 2 uncertainties. The sample of 
protostars we have analyzed has an estimated maximum 
mass m„ = SM©; a larger sample would presumably in- 
clude more massive stars and enable a stronger test of the 
theories. If the sample were sufficiently large, one would 



be able to directly determine the role of large, FU Ori 
type outbursts on the growth of protostars. Study of the 
very faint protostars, the VELLOs, should determine the 
relative proportion of those that are in a quiescent state 
between outbursts and those that are faint because they 
have very low mass (as we have assumed). A monitoring 
program would permit one to characterize the variability 
of the protostars, which we have simply taken as increas- 
ing the dispersion in the PLF by a factor of two. Finally, 
more accurate measurements of the physical conditions 
in the star-forming clumps would enable more accurate 
theoretical predictions of the accretion histories of the 
protostars in the sample. 
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APPENDIX 

PLF ALTERNATIVE DERIVATION 

The PLF may be obtained either by integrating '^p^iL^m) over m/ or over m, where in either case the resulting 
PLFs are identical. Above we exclusively use the former formalism. For completeness, we give the latter PLF definition 
here: 



dlnm'^ p2{L , m), 

ipp2 [m,mf{L,m)] 



dlnm ■ 



dhiL 



(Al) 
(A2) 



where the lower limit of integration is given by equation (5) with m = m(L, m/). Note that this formulation does not 
work when the luminosity is independent of the final mass, as in the case of Isothermal Accretion. 

THE PLF FOR CONSTANT RADIUS 

Accurate evaluation of the PLF requires allowing for the dependence of the protostellar radius on the mass and 
accretion rate. However, the radius is almost always within a factor 2 of r = 2.5 Rq, and if we take r to be constant 
the analysis is simplified considerably. 

Combining equations (1) and (20), we express the accretion luminosity as 

r 1/2 
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is the luminosity for r = 1 Rq and m = m/ = 1 Mq in the untapered case. If we let 

L{1) 

(where r is in units of Rq), then the relation of m and m/ to L is 
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Isothermal Sphere Accretion 
In this case, we have m = ^, so that equation (29) gives 

*p(L) = Vp(m = = -/ dlnmf^imf). (B6) 

Tapered Isotherm.al Sphere Accretion 
In this case, equation (B5) for the accretion luminosity gives 

ruf 

which can be solved for mf{£, m), 

mf{i,m) = (B8) 

Note that m> £. The result for the PLF is then (eq. A2), 



1 - (i/my 



(B9) 



The minimum possible value of m for a given luminosity corresponds to w/ = As m increases, m/ decreases until 
it reaches a minimum, 

/33/2\ 

min(m/) =( — ]£, (BIO) 



and it then increases again; as a result, the maximum possible value of m for a given luminosity also corresponds to 
rrif = m„. The limits of integration of the PLF, mmin and rumax, are therefore given by the roots of equation (B7) 
with m/ = m„ and that satisfy the condition m > £. In order for mmin and rumax to be less than m„, it is necessary 
that min(m/) be less than m„, which sets an upper limit on the luminosity, 

4 = (,^) m„. (Bll) 



33/ 

One can show that iu is the maximum luminosity allowed by equation (B7) with my = m„. It is also necessary that 
rUf exceed me, which becomes an issue when £ < 2me/3^/^, so that min(m/) < m^. In this case, the range of mass 
between the roots of equation (B7) with ruf = mg that satisfy m > £ is excluded. The resulting range of integration 
extends from the smaller root of equation (B7) with m/ = m„ to the smaller root of this equation with m/ = me, and 
then from the larger root equation (B7) with mf = me to the larger root of the same equation with mf = m„. 

Untapered Turbulent Core and Competitive Accretion 
Since for untapered Turbulent Core and Competitive Accretion, the luminosity equation (B5) implies 

1/0/ -i) / £ \ 



it follows that 
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(B13) 



The minimum mass for a given luminosity corresponds to the solution of equation (B12) with m/ = m„. 



mi'-' 



mmin = ——■ . (B14) 



There are two conditions that set the maximum value of m: First, the minimum value of m/ is m^, so that 



(B15) 



Second, the requirement that m be no more than the final mass, m/, implies 

m^^<£^/^^+ifK (B16) 

The correct value of mmax is the lesser of these two values. The upper limit on the luminosity occurs when m = m/ = 

rriu, 

£u = mi+''. (B17) 

The condition £ < £u ensures that mmax < m„. 
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In this case, we have 



which leads to 
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Tapered Competitive Accretion 
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Evaluating the partial derivative dlnL/dlnnif, we find that the luminosity function is 
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Just as in the case of tapered Isothermal Sphere accretion, there is a minimum value of m/ as a function of m for a 
given luminosity, 

•1111/4^ 
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This is less than provided the normalized luminosity is less than 
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For £ < £u, the limits of integration, TOmin and mmax. are the roots of equation (B18) with m/ = m„ and < m < rn„. 

Tapered Turbulent Core Accretion 



In this case, equation (B5) becomes 



so that 

The luminosity function is then 
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The limits of integration are the solutions of equation (B23) with mf = m^ and < m < m„, just as in the case 
of tapered Competitive Accretion. The maximum possible luminosity can be found directly from maximizing £ in 
equation (B23), or, equivalently, by evaluating min(m/) and requiring that it be less than m„: 
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Two-Component Turbulent Core Accretion 
For Two-Component Turbulent Core Accretion, the accretion rate is 

m = rhis ^1 + T^^mm/i/^^ 

In this case, we define the normalized luminosity as 

rL / / \ 1/2 

e = j;^ =m[l + nlmmf'/^j , 

so that the final mass for a protostar of mass m and normalized luminosity £ is 
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Evaluating 9 In L/c? In m/ , we find that the PLF is 



.2 \ 2' 
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(B30) 



The lower limit of integration, mminj is the solution of equation (B27) with nif = to„. Unless the luminosity is very 
low, the upper limit of integration is set by the condition m = m/, so that mmax is the root of the equation 

f = mn^J" (l + 7e^mma^3/2) . (B31) 

However, m/ cannot be less than the minimum stellar mass me, so that in general mmax is the solution of the equation 



'2 _ 2 
m-max 



1 + 7^^m„laxmax(mmax^^^,my^) . (B32) 
The maximum possible luminosity occurs when m = ruf = m„, so that 

£u=ml(l+Tllml/^). (B33) 

For £<iu, both 

Tlinin and TTlmax less than m„. 

Tapered Two- Component Turbulent Core Accretion 

In terms of 

g=(l + nlmmf'/^y^\ (B34) 

the equation for tapered accretion is 

m = rhisg (i - • 
Integrating this equation gives the relation for the age as a function of mass, 

t-& = ^^ 172(5-1)- 

2tf mIs7^^m/V2 

This relation shows that the formation time with tapering, tf, is twice the value without, tf = 2tfo with 

2 

mIs7^^m/ 

The tapering factor is then 



*/o= ^..^4^.V2 (5/-l)- (B37) 



so that the normalized luminosity is 



£ = mg (j^^ ^ ■ (B39) 

The protostellar luminosity function is then given by equation (29) with m/ determined numerically from equation 
(B39) and with 

Based on comparison with the cases of the tapered IS and tapered TC cases, the upper and lower limits of integration 
for the PLF are found from numerical solution of equation (B39) with ruf = m„. 
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